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Preface 


Finally,  as  all  our  observations,  on  account  of  the  imperfection  of  the 
instruments  and  of  the  senses,  are  only  approximations  to  the  truth, 
an  orbit  based  only  on  the  six  absolutely  necessary  data  may  be  still 
liable  to  considerable  errors.  In  order  to  diminish  these  as  much  as 
possible,  and  thus  to  reach  the  greatest  precision  attainable,  no 
other  method  will  be  given  except  to  accumulate  the  greatest  number 
of  the  most  perfect  observations,  and  to  adjust  the  elements,  not  so 
as  to  satisfy  this  or  that  set  of  observations  with  absolute  exactness, 
but  so  as  to  agree  with  all  in  the  best  possible  manner.  [1 :51] 

Carl  Friedrich  Gauss,  quoted  above,  was  the  founder  of  estimation  theory.  He  laid 

the  ground  work  for  all  modern  estimation.  Today,  optimal  estimation  has  been 

claimed  by  the  Kalman  filter.  This  study  used  the  iterated,  extended  Kalman  filter 

to  investigate  the  performance  of  an  on-board  filter  used  to  estimate  the  relative 

positions  of  a  cluster  of  satellites  navigating  in  the  presence  of  real  dynamics.  But 

it  wasn’t  just  the  Kalman  filter  that  allowed  this  research  to  be  possible. 

Thank  you  to  all  who  aided  me  in  this  endeavor.  I  wish  to  specifically  thank 

Dr.  William  E.  Wiesel,  my  thesis  advisor,  for  his  patience,  guidance  and  good 

humor  during  this  research.  Additionally,  I  want  to  thank  Dr.  Wiesel  and  Dr. 

Rodney  Bain  for  the  enjoyment  I  received  from  their  astrodynamics  classes.  Also 

I  thank  Dr.  Peter  Maybeck  for  his  time  and  counsel.  I  wish  to  thank  Capt  John 

Gustafson  for  his  advice. 
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My  unending  love  and  thanks  are  due  my  wife  Joan  and  son  Stephen  for 
their  patience,  understanding  and  support  for  these  past  months.  Finally,  1  thank 
our  Lord  for  the  endurance  to  see  this  through.  I  pray  in  some  way  this  serves  His 
purpose,  otherwise  it  will  have  been  for  nothing. 


J.  Timothy  Middendorf 
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Abstract 


Previous  work  in  the  area  of  estimation  of  relative  positions  within  a  satellite 
cluster  showed  favorable  results.  However,  the  work  was  done  using  point  mass 
orbits  in  the  t.uth  model.  This  thesis  investigates  the  estimation  of  relative  satellite 
positions  operating  in  near  circular  orbit  including  the  Jj  term  in  Earth's 
geopotential.  The  iterated,  extended  Kalman  filter  is  used  as  the  on-board 
estimator  in  order  to  gain  better  performance  in  the  face  of  the  non-linearities.  The 
dynamics  in  the  estimator  are  based  on  the  Clohessy-Wiltshire  equations  for 
relative  orbital  motion.  Inputs  to  the  host  estimator  are  range  measurements  from 
each  satellite  in  the  cluster.  The  relative  position  to  the  host  satellite  is 
investigated.  A  comparison  of  the  true  error  and  the  rms  covariance  was 
performed. 
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NAVIGATION  OF  A  SATELLITE  CLUSTER 


WITH  REALISTIC  DYNAMICS 


I.  Introduction 

Past  studies  have  shown  the  feasibility  of  using  a  recursive  filter  for 
on-board  estimators  to  determine  the  relative  positions  of  satellites  for  use  as  a 
space  based  radar  system  [5,7,10].  Good  results  were  obtained  using  two-body 
point  mass  astrodynamics  and  an  on-board  estimation  model  with  the  dynamics 
based  on  the  Clohessy-Wiltshire  equations. 

The  initial  work  at  AFIT  was  done  by  Captain  Michael  Ward  using  the  U-D 
covariance  factorization  Kalman  filter  with  the  Clohessy-Wiltshire  equations. 
However,  he  encountered  observability  problems  in  the  down  range  state 
component  [10:2-17].  Further  investigation  into  the  observability  problem  was 
conducted  by  Captain  Sherrie  Norton  Filer  [5]  and  Captain  Stephen  Johnston  [7], 
Captain  Filer  obtained  mixed  results  from  her  investigation  into  unobservable 
states.  Captain  Johnston  furthered  the  effort  and  realized  the  Clohessy-Wiltshire 
equations  could  not  determine  the  position  and  velocity  of  each  satellite  from  the 


1 


origin  on  the  reference  orbit  based  on  range  measurements  between  the  satellites. 
Only  relative  position  and  velocity  could  be  determined  from  the  host  satellite  to 
each  other  satellite  [7:29].  The  host  satellite  considers  itself  the  center  of  the 
coordinate  frame.  After  this  was  identified  and  analysis  done,  the  results  were 
favorable.  He  obtained  3  cm  accuracy,  well  beyond  the  25  m  requirement. 

However,  the  truth  dynamics  were  based  on  a  homogeneous,  spherical 
Earth;  effectively  treating  the  Earth  as  having  point  mass  gravity.  This  ignores  one 
of  the  most  dominant  perturbative  influences  of  a  satellite’s  orbit.  The  Earth’s 
rotation  produces  a  bulge  at  the  equator  and  large  structures  such  as  mountain 
ranges  produce  variations  in  the  local  mass.  The  effect  of  this  oblateness  and 
non-homogeneity  is  a  non-uniform  gravitational  field.  Although  other  effects  are 
present  as  well,  such  as  air  drag  and  solar  pressure,  this  thesis  employs  a  truth 
model  which  accounts  for  the  non-uniform  gravity  field  only. 

This  thesis  investigates  the  estimation  of  relative  satellite  motion  in  the 
presence  of  these  real  dynamics.  The  truth  model  computer  program  was 
developed  to  include  the  J2  term  of  Earth’s  geopotential.  Since  there  are 
significant  non-linearities  between  the  fruf/?  dynamics  and  the  estimator  dynamics, 
the  estimator  developed  was  the  iterated,  extended  Kalman  filter  [9:58], 

The  nominal  orbital  altitude  (1000  km),  cluster  radius  (500  meters),  and 
accuracy  requirements  (25  meters)  remain  the  same  as  the  previous  investigations 
[10:1-1  ,App  A).  The  accuracy  is  based  on  the  requirement  of  the  radar  to  form  a 
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cohesive  image  using  the  combined  data  from  each  satellite.  In  order  to  achieve 
this,  the  relative  position  of  each  satellite  must  be  known  to  at  least  one  quarter 
of  the  radar  wavelength  [10:1-1]. 

The  initial  concept  employed  a  cluster  of  up  to  10  satellites  orbiting  in  near 
circular,  low  Earth  orbit.  Capt  Johnston  realized  each  satellite’s  relative  state  was 
independent  of  the  other  [7:29].  Therefore,  only  a  two  satellite  cluster  will  be 
investigated.  The  procedure  can  be  generalized  to  include  more  satellites  in  the 
cluster. 

In  this  work,  the  origin  of  the  rotating  reference  frame  is  centered  on 
satellite  #  1  which  will  contain  the  on-board  estimator  being  investigated.  The  truth 
model  is  the  ruler  against  which  the  estimator  is  placed.  A  comparison  of  true 
error  and  rms  covariance  will  be  made  [8:337]. 
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II.  Background 


The  fundamental  purpose  of  this  investigation  is  to  compare  the  estimation 
with  real  dynamics.  This  is  done  by  generating  a  truth  model  based  on  real 
dynamics  and  comparing  with  the  results  achieved  from  the  filter’s  estimation  of 
the  state.  In  this  way,  the  same  mission  can  be  run  over  and  over  with  different 
random  noise  in  the  range  to  analyze  the  filter’s  performance. 

The  truth  model  computer  program  was  written  to  include  the  Jj  term  of 
Earth's  geopotential.  The  J2  term  is  the  most  dominant  term  in  the  geopotential, 
beyond  the  Newtonian  point  mass  orbit  potential.  It  accounts  for  both  periodic  and 
secular  variations  in  an  orbit.  The  Jg  effects  cause  the  secular  variations  to  be  the 
dominant  perturbation  source  in  the  Mean  Anomaly,  M,  Argument  of  Perigee,  co, 
and  the  Right  Ascension  of  the  Ascending  Node,  fl,  near  1000  km.  Figure  1  and 
Figure  2  show  the  Newtonian  (constant)  value,  the  theoretically  predicted  secular 
variation  and  the  actual  secular  variation  in  co  and  Q  respectively. 

The  predicted  secular  variations  are  given  by  [4:369]: 


0)  =  (Oq 

Q  =  Q, 


(b?  =  COg  + 


2a2(i  - 


3T|J2R|f 
23^(1  -  e")' 


cost 


(1) 
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Effect  of  Oblateness  on  the  Argument  of  Perigee 


I _ Time  ITUs] _ 

Figure  1 .  Comparison  of  the  Newtonian,  predicted  and  actual  secular  variation 
in  the  Argument  of  Perigee,  co,  due  to  Jj. 


Figure  2.  Comparison  of  Newtonian,  predicted,  and  actual  secular  variation  in 
the  Right  Ascension  of  the  Ascending  Node,  Q,  due  to  Jj. 


5 


where: 


Ti  =  Mean  Motion 

Jg  =  Zonal  harmonic  coefficient  of  the  geopotential 
a  =  Semi-major  axis 
e  =  Eccentricity 
I  =  Inclination  of  the  orbit 

The  actual  secular  variation  shows  oscillations  periodic  with  the  orbit.  This 
is  due  to  the  periodic  variations  in  a,  e  and  t  of  the  orbit  due  to  J2.  These  plots 
verify  the  proper  functioning  of  the  orbit  determination  part  of  the  truth  model.  The 
actual  variations  with  the  periodicity  in  the  orbit  follow  the  predicted  variations  as 
expected. 

The  estimator  computer  program  was  written  as  an  iterated,  extended 
Kalman  filter  employing  the  Clohessy-Wiltshire  equations  for  the  dynamics  to 
predict  the  state.  It  uses  range  measurements  from  the  truth  model  to  update  the 
predicted  state.  The  truth  model  generates  the  relative  position  and  ve'ocity  of  the 
cluster  and  provides  the  estimator  with  inputs  of  range  corrupted  with  white 
Gaussian  noise.  The  estimator  uses  these  inputs  to  generate  its  correction  to  the 
relative  position  elements  of  the  state.  It  also  generates  an  estimate  of  its  error, 
the  covariance.  This  covariance  is  compared  to  the  true  error  to  allow  tuning  the 
filter. 
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2. 1  Truth  Model 


2.1.1  Coordinate  Frame  Definition 

The  problem  is  presented  using  three  coordinate  frames.  Two  of  these 
frames  will  be  employed  through  out  the  thesis  and  one  serves  only  to  ease  the 
set-up  of  the  initial  conditions  of  the  second  satellite.  The  frames  are;  the  inertial 
frame  centered  on  the  Earth,  the  rotating  frame  centered  on  the  host  satellite,  and 
the  third  is  the  inertial  frame  aligned  with  the  rotating  frame. 

The  first  frame  is  referenced  in  the  truth  model.  The  equations  of  motion 
in  the  truth  model  for  each  satellite  are  defined  in  the  inertial  frame  with 
components  (  X,  Y,  Z  ).  X  is  aligned  with  the  First  Point  of  Aries  (T),  Z  the 
rotation  axis  of  the  Earth,  and  Y  completes  the  right  hand  orthogonal  set. 

The  second  frame,  the  rotating  frame,  is  fixed  to  satellite  1.  It  has 
components,  (  f,  i,  h  )  where,  r  is  along  the  radius  vector  from  the  center  of 
the  Earth  to  satellite  1,  h  is  perpendicular  to  the  orbital  plane  and  i 
completes  the  right  hand  orthogonal  set.  This  frame  is  the  frame  referenced  in  the 
estimator.  Once  satellite  2’s  position  is  found  inertially,  the  relative  position  from 
satellite  1  is  found  and  transformed  into  the  rotating  frame. 

The  third  frame  is  only  used  to  define  the  initial  inertial  conditions  of  satellite 
2  using  the  Clohessy-Wiltshire  equations  to  provide  appropriate  initial  velocities. 
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This  frame  is  simply  the  inertial  frame  (  X,  Y,  Z  )  rotated  to  align  with  the  rotating 
frame  and  has  components  {  P,  Q,  W  ).  After  the  initial  conditions  are  specified 
only  the  inertial  (  X,  Y,  Z  )  and  the  rotating  (  f,  i,  h  )  frames  will  be  used. 
Figure  3  illustrates  these  frames. 


Figure  3.  Coordinate  frame  definition. 
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2.1.2  Initial  Conditions 


Initial  conditions  are  specified  by  a  reference  orbit  for  the  host  satellite  and 
specifying  the  relative  state  of  satellite  2.  Once  a  reference  orbit  is  specified  for 
satellite  1 ,  the  host,  satellite  2’s  initial  relative  state  is  defined  in  the  rotating  frame. 
Satellite  1  ’s  initial  state  is  defined  in  the  inertial  frame  and  expressed  in  the  rotating 
frame.  This  initializes  satellite  1’s  state  in  the  truth  model  and  the  filter.  Satellite 
2’s  initial  state  is  defined  in  the  rotating  frame  about  satellite  1 .  Its  state  is  also 
expressed  in  the  inertial  frame.  This  initializes  both  the  truth  and  filter  models  for 
satellite  2. 

The  initial  position  and  velocity  of  satellite  1,  and  therefore  the  rotating 
reference  frame,  are  arbitrarily  chosen  to  have  the  following  classical  elements, 
a=1. 156784906  (1000  km  altitude),  e=0.00011, 1=28.5°,  q)=30°.  Q=50°  and  v=0°. 
The  initial  location  of  this  rotating  frame  in  the  inertial  frame  (  X,  Y,  Z  )  is  easily 
obtained  from  the  classical  elements  [1 :71 ,82]  via  the  (  P,  Q,  W  )  frame. 

The  position  and  velocity  in  the  P,  Q,  W  frame  are  expressed  as; 


r  =  rcosvP  +  rsinvQ 


V 


f  P 

,a(1-e2) 


\1/2 

[-sinv 


P 


+  (e+cosv)Q] 


(2) 

(3) 
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These  are  rotated  into  the  inertial  frame  {  X,  Y,  Z  )  by: 


rx 

ry 

=  R 

” 

and 

r 

X  > 
.>  > 

=  R 

Vp 

Vq 

Vz 

Vw. 

(4) 


where, 


R  = 


cosQcosco-sinQsinfocosl 

sinQcosco+cosQsinfocost 

sincDSinl 


-cosQsinco-sinfl  COSO)  cost 
-sinQsino)+cosQcos(ocosl 
coso)sinl 


sinQsinl 

-cosQsinl 

cost 


(5) 


With  the  classical  elements  given  above,  the  initial  position  and  velocity  of 
the  rotating  frame  expressed  in  the  inertial  frame  in  canonical  units  are: 


Ri 


0.25453861 

1.09403658- 

0.27595466 


DU 


(6) 


V, 


-0.84098533 
-  0.09874202- 
0.38425099 


DU 

TU 


(7) 


where,  1  DU  =  6378.145  km  and  1  TU  =  806.81 18744  secs.  These  are  found  by 
defining  Earth's  gravitational  constant,  =  1  DU^/TU^  =  3.986012  X  10^ 

km^/sec^. 
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Since  the  rotating  frame  is  in  a  near-circular  orbit  (e=. 00011),  the 
transformation,  [t]  ,  from  the  inertial  frame  (  X,  Y,  Z  )  to  the  rotating  frame 
(  f,  I,  h  )  is  determined  from  the  host  satellite’s  position  and  velocity: 


'  N 

f 

X 

T, 

r 

X 

i  =[t1- 

h 

Y 

Z 

'  — 

T* 

Tn 

Y 

Z 

(8) 


where, 


[T.l 

[TJ 

(Tn) 


RiJ 


[x.Y,2] 


=  (TJ  X  [T,] 


R  X  V 


L|r  X  v|J[jc>,2] 


(9) 


and  R  and  V  are  the  radius  and  velocity  of  satellite  1  in  the  inertial  coordinates. 

The  position  elements  of  satellite  2’s  state  are  randomly  placed  about  the 
rotating  origin: 

^2  =  (500  meters)  (10) 


where  is  a  vector  of  random  numbers  between  -0.5  and  0.5. 
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It  is  rotated  into  the  inertial  frame  (  X,  Y,  Z  )  to  be  used  by  the  truth  model.  The 
position  of  satellite  2  in  the  inertial  frame  is  found  by: 


R 


‘2[x.y,2] 


MtI 


(11) 


With  the  position  of  satellite  2  specified  in  both  frames,  the  velocity  is 
defined  using  the  Clohessy-Wiltshire  equations.  The  Clohessy-Wiltshire  equations 
(see  Section  2.2.1)  allow  appropriate  values  for  radial  and  out-of-plane  velocity 
components  in  the  (  P,  Q,  W  )  frame.  Both  the  radial  and  out-of-plane  velocity 
components  are  found  by  multiplying  the  relative  distance  expression  by  the  mean 
motion  of  the  reference  orbit  and  then  equating  to  the  relative  velocity  expression. 
Since  this  frame  is  initially  aligned  with  the  rotating  reference  frame, 

VjP^Tir^P  (12) 

V2W=Tir2W  (13) 


where  ri  =  the  mean  motion  of  the  reference  trajectory. 

In  order  to  maintain  cluster  integrity  by  reducing  drift,  the  orbital  periods  of 
each  satellite  must  be  the  same.  The  orbital  period  in  the  two-body  Newtonian 
point  mass  orbit  is  a  function  of  the  semi-major  axis  only  and  completely  removes 
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drift.  Letting  this  be  the  starting  point  for  the  perturbed  orbit,  the  two-body  energy 
equation  is  used  to  determine  the  third  component  of  the  velocity.  Utilizing  this 
equation  and  assuming  the  semi-major  axis  is  constant  and  the  geopotential  is 
based  on  a  point  mass  yields: 


f  \ 

V2  Q  =  ■ 

f' 

1  _  1 

-(V.-V,)P  -(V.-V,)W- 

These  velocity  components  are  expressed  in  the  rotating  frame  to  be  used 
by  the  estimator  as  follows: 


®  ~  (O)  X  r2)[f ; 


(15) 


where,  initially.  =  (^2  "  V,  )[p.,<5,yA/]  • 


The  velocity  components  are  expressed  in  the  inertial  frame  to  initialize  the  truth 
model  by 


^2[x.Y.Z]  ~  .1"  ]  ^2[f>.6,A] 


(16) 
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Therefore,  the  initial  state  to  the  estimator  is 


x,(0) 


and  the  initial  state  to  the  truth  nodel  is 


Xt(0) 


R, 

V, 

R, 

Va 


[x.y,2] 


(17) 


(1b) 


2.1.3  Dynamics 

Having  defined  the  initial  conditions,  the  truth  model  state  must  be  found  at 
subsequent  times.  Previous  studies  used  Newtonian  point  mass  dynamics  for  the 
truth  model.  As  addressed  above,  this  maintained  the  cluster’s  integrity  well. 
However,  the  Earth  is  non-spnerical  and  non-homogeneous,  which  causes  drift  to 
appear  in  the  real  dynamics.  Orbits  based  on  the  point  mass  relative  dynamics 
and  the  Clohessy-Wiltshire  dynamics  retrace  themselves  since  there  is  no 
perturbations  present.  Figure  4  is  a  plot  of  10  orbits  of  both  the  point  mass 
relative  dynamics  and  the  Clohessy-Wiltshire  dynamics.  There  is  no  perceptible 
difference  after  10  orbits.  Figure  5  shows  the  distinct  difference  in  the  relative 


14 


Z  [meters]  Z  [meters]  X  [meters] 


500 


TRUTH  and  CLQHESSY-WILTSHIRE  Relative  Orbits  without  J2 


-500  400  -300  -200  -100  0  100  200  300 


X  (meters  1 


Figure  4.  Point  mass  and  Clohessy-Wiltshire  orbits  (10  orbits). 
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Z  [meters]  Z  [meters]  X  [meters] 


Y  [meters] 


Figure  5.  flea/ dynamics  and  Clohessy-Wiltshire  orbits  (10  orbits). 
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orbits  of  the  real  dynamics  and  the  Clohessy-Wiltshire  dynamics.  Later  this  will 
play  an  important  factor  in  the  performance  analysis. 

This  difference  is  due  to  the  oblateness  of  the  Earth  and  the  gravitational 
potential.  The  Earth’s  geopotential  including  the  zonal  harmonic  term  and 
neglecting  the  other  non-vanishing  higher  order  terms  is  [12:84]: 

r 


U  =  ii  -  ii 
R  R 


R© 


l(3cos2e  -  l)jj 


(19) 


where: 


COS0  =  — 

R 

Rq  =  radius  of  the  Earth 
[L  =  Earth's  gravitational  constant 


I 


(20) 


The  truth  model  uses  the  potential  of  Eq  (19)  expressed  in  the  inertial  frame  as 


p.  _ 

f 

3Z2 

X^+Y^+Z^ 

\ 

-  1 

y 


(21) 


With  set  to  0,  the  geopotential  is  simply  the  Newtonian  point  mass  potential. 
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The  future  position  and  velocity  of  satellites  1  and  2  are  determined  inertially 
by  numerically  integrating  the  equations  of  motion.  Integration  is  done  using  a 
predictor-corrector  routine  called  NAMING  [12:110]  (see  Appendix  A).  NAMING 
requires  a  subroutine  called  RNS  (Right  Nand  Side).  This  subroutine  contains  the 
right  hand  side  of  the  state-differential  equations  of  motion: 


X  =  V, 
Y=V, 
Z  =  V, 


X  =  V, 
Y=V, 
Z  =  V, 


_  au  _ 

_pX 

1 

3  ^2^® 

5Z2 

Y 

-  1 

ax 

R"* 

2  R2 

[r^ 

)- 

- 

( 

Y 

_  au  _ 

_pY 

1 

_  3  J2R4 

5Z== 

- 1 

■  aY  ■ 

R=* 

2  R2 

[r^ 

). 

- 

f 

Y 

-  au  _ 

_ 

1 

3  ^2^ 

3  - 

5Z2 

‘  az 

R" 

2  R2 

R2 

(22) 


A  sample  rate  for  the  estimator  is  preselected  in  order  to  generate  range 
data  for  the  estimator.  The  sample  rate  determines  how  frequently  the  states  of 
satellites  1  and  2  will  be  used  to  calculate  the  range  provided  to  the  estimator.  At 
each  of  these  sample  times  the  relative  state  for  satellite  2  is  determined  by: 


^r8l[x.?,2]  ~  {^^2(0  ^1(0} 

>  {V3(()  -  V,(o} 
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The  coordinate  transformation  matrix,  [t]  ,  is  re-evaluated  by  Eq  (9).  The 


relative  position  and  velocity  are  rotated  into  the  rotating  frame  as  follows; 


where, 


r 

V  -  (co  X  7} 


V 


rel[X,'i',2] 


(24) 


(25) 


Finally,  the  truth  model  generates  the  range  measurement  for  the  estimator 
by 

Zn  =  ^  (26) 

where  u^  is  a  zero-mean,  white  Gaussian  noise  with  its  associated  covariance  of 
Rn  [8:330].  The  noise  is  the  representation  of  measurement  inaccuracies  as  well 
as  other  error  sources  [10:2-6]. 

2.2  Filter  Model 

The  choice  of  the  iterated,  extended  Kalman  filter  was  made  in  light  of  the 
significance  of  the  non-linearities  of  the  system,  both  in  the  observation  relationship 
and  the  dynamics  [9:58].  The  iteration  occurs  at  each  data  sample  in  order  to 
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modify  the  reference  trajectory  with  a  better  estimate  of  the  state.  After 
convergence,  which  was  chosen  to  be  a  difference  of  0.0001  in  each  iterated 
position  element,  or  a  specified  number  of  iterations,  chosen  to  be  8,  the  final 
reference  state  is  deemed  the  new  estimate  with  its  associated  covariance. 

One  poiential  problem  with  the  Kalman  filter  is  "starting"  the  filter,  since  the 
first  time  propagation  occurs  before  any  data  is  sampled.  Therefore,  the  initial 
state  is  considered  known  from  the  truth  model  as  described  in  Section  2.1.2.  The 
initial  covariance,  Pq,  must  also  be  specified.  Pq  will  be  considered  diagonal  with 
initial  position  and  velocity  covariance  of  25  rrf  and  1  (m/sf,  respectively.  These 
result  from  an  assumption  of  an  initial  position  accuracy  of  5  m  and  velocity 
accuracy  of  1  m/s.  These  will  be  replaced  with  more  appropriate  values  as  the 
filter  is  tuned. 

2.2. 1  Dynamics  (Clohessy-Wiltshire  Equations) 

The  state  is  propagated  forward  to  the  time  of  next  data  acquisition.  Since 
the  dynamics,  governed  by  the  Clohessy-Wiltshire  equations,  are  linear  (see 
below),  the  state  and  covariance  propagation  are  defined  by  [8;220]; 

Ht:)  =  MC.)  (2^) 
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P(^-)  =  -  G.iUQ^iUGjiU 


(28) 


where, 

Qjj  is  the  covariance  of  the  dynamics  noise. 

x(/,')  is  the  estimated  state  after  time  propagation. 

x{t-)  is  the  estimated  state  after  the  measurement  update. 

Gd  is  the  identity  matrix  for  this  equivalent-discrete-time  representation  of 
a  continuous-time  system  [8:377]. 
cf)  is  the  state  transition  matrix. 


The  state  is  propagated  to  time  (ti)  by  means  of  O: 


<D  = 


ax(0 

ax(fo) 


(29) 


Figure  6  illustrates  the  state  propagation  and  data  acquisition  [8:207]. 
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The  equations  of  relative  motion  used  to  develop  the  Clohessy-Wiltshire 


equations  are  defined  as  follows  [13:80]: 


5f  -  2riro59  -  3ri^6r  =  0 


(30) 


ro6e  +  2Ti5r  =  0 


(31) 


5z  +  ti^5z  =  0 


(32) 


where,  6r,  ro50,  and  5z  are  the  f,  i,  h  coordinates  respectively. 


These  coordinates  will  be  defined  as  x,  y,  and  z  for  the  remainder  of  the  thesis. 
It  should  be  noted,  these  are  a  first  order  form  of  the  Newtonian  point  mass 
relative  equations  of  motion. 

Eqs  (30),  (31),  and  (32)  are  readily  integrated  since  they  are  a  set  of  linear, 
constant-coefficient  differential  equations.  Applying  the  initial  conditions. 


^0 

^0 

x(o)  =  • 

Vo 

x(o)  =  ' 

Vo 

,^o  , 

^0  . 
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the  solutions  are  [13:80-81]: 


z(/)  =  ZqCOS  r]t  +  iisin  r\t  (35) 

n 

x{0  =  (2yo  +  STixJsin  r\t  +  XpCos  r\t  (36) 

y(0  =  (-3yo  -  6tiXo)  +  (6tiXo  +  4y(,)cos  T]t  -  Zx^sln  T]t  (37) 

z(t)  =  -ZoTisin  rif  +  z^cos  rif  (38) 
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Now  the  state  transition  matrix  is  determined  by  Eq  (29): 


4-3cos\y 

0 

0 

Isinvy 

1(1  -cosvy)  0 

6(sin\|r-v) 

1 

0 

i(cosv-l) 

lsin\j/-i.y  0 

0 

0 

COS\|/ 

0 

0  -Lsinw 

n 

3Tisin\j/ 

0 

0 

cosy 

2sin\|/  0 

6Ti(cos\i;-1 ) 

0 

0 

-2sin\|/ 

-3+4cos\y  0 

0 

0 

-Tisinxy 

0 

0  cosy 

where,  v  =  t]/  and  /  is  the  sample  time. 

The  state  is  propagated  to  the  new  data  acquisition  time  by  Eq  (27),  and  in 
general,  Eq  (27)  becomes: 


cD^  0  -  0 

MCr) 

x(f;)  = 

o 

...  e 

u 

...  o 

r 

O 

O 

e 

m 

i 

i _ 

(40) 


where  the  subscripts  reference  each  satellite  in  the  cluster. 
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<1)  need  only  be  evaluated  once  given  the  sample  interval,  (f,.,)  to  (t,)  .  This  is 
one  of  the  major  benefits  of  linear  dynamics. 

2.2.2  Iterated,  Extended  Kalman  Filter 

With  the  initial  state  and  covariance  propagated  to  time  (//),  the 
measurement  update,  is  incorporated  by  means  of  the  state  and  covariance 
update  equations.  Data  is  typically  not  a  scalar,  and  therefore  the  expectation  of 
the  data  is  typically  a  vector  and  is  denoted  by  h.  The  bold  notation  is  retained 
for  both  the  data  and  its  expectation  for  the  remainder  of  the  thesis.  Since  z,  is  the 
range,  it  is  expected  to  be  of  the  form  of  h,  which  is  defined  as 

h  =  {xl  ^  y|  +  zlY'^  (41) 


In  general,  h  becomes 


' 

;  > 

.{x.' 

*  Vs 

•  4}"’ 

‘’.I 

(42) 


where,  x,  y,  and  z  are  the  components  in  the  (  f,  i,  h  )  ir-ame. 
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The  expectation,  h,  is  linearized  and  evaluated  at  x(ti)  by 


(43) 


which  becomes 


H  = 


X 

h 


y  £  0  0  0 

h  h 


(44) 


In  general,  H  has  the  form 


0  0  -  0 
0  H,  0  •••  0 

:  :  :  0 

0  0  0  -  H,., 


(45) 


where. 


H 


S-1 


Z 


s-1 


h 


t-i 


0  0  0 


(46) 
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The  state  update  is  given  by  [9:59]: 


Xk.i  =  x(f')  +  K(t,){z,.  -  h(x„f^)  -  H(x,.f.)[x(/,')  -  X,]} 


(47) 


where, 

K((,)  =  P((-)H'(x„/,)[H(x,.t,)P(f-)H^V/,)  *  R(0]'’ 

fork  =  0,  1,  2 . N-1. 

R(f,)  is  the  covariance  of  the  noise  in  the  data. 

The  iterated,  extended  Kalman  filter  update  equations  are  a  slightly  modified 
form  of  the  extended  Kalman  filter  update  equations  [9:44,59].  The  iterated 
method  uses  the  updated  state  estimate  x{f/')  ,  generated  by  the  standard 
extended  Kalman  filter,  "as  a  better  state  estimate  than  x  ( f,' )  for  evaluating  h  and 
h  in  the  measurement  update  relations"  [9:58]. 

This  is  important  since  the  h  vector  is  the  expectation  of  the  data,  z,,  and 
is  used  to  evaluate  the  residual. 


r,  =  z,  -  h(x(0) 


(49) 
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found  in  the  state  update  equation,  Eq  (47). 


It  is  seen  that  x,  is  actually  x{tj)  given  by  the  standard  Kalman  filter 
[9:59].  When  the  specified  nurfiber  of  iterations  have  been  made  or 
when  I  -  x^.,  |  is  less  than  a  predetermined  amount,  x^  is  declared  the  new 
reference  trajectory:  that  is, 

xit!)  =  x^  (^0) 


The  new  reference  trajectory  is  now  permanently  updated,  with  a  covariance 
of 


P(/;)  =  P(/-)  -  K(gH  x(f,'),dP(f;) 


(51) 


This  iterative  method  reevaluates  the  observation  relations,  h  and  H,  with 
each  x^  to  achieve  a  better  reference  trajectory.  This  done  to  cc'rect  for  the 
nonlinearities  in  the  observation  and  the  dynamics  of  the  system. 
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III.  Performance  Analysis 


3. 1  Performance  Evaluation 

Assessment  of  the  filter’s  performance  is  accomplished  by  comparing  the 
true  error  and  root  mean  square  covariance  [8:337].  By  graphically  comparing  the 
true  error  and  the  root  mean  square  covariance,  the  filter  can  be  "tuned".  This 
method  of  tuning  will  give  a  general  assessment  of  the  filter’s  capabilities  of 
tracking  the  other  satellites  in  the  cluster. 

In  order  to  perform  the  comparison,  the  position  elements  of  both  models 
are  extracted  from  the  complete  state  by 


y,  =  c.x, 

y  =  Cx 


(52) 


Since  the  states  are  the  same  for  both  models. 


C  =  C, 


1  0  0  0  0  0 
0  1  0  0  0  0 
0  0  1  0  0  0 


(53) 
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The  true  error  is  depicted  in  Figure  7  and  is  defined  as 


e,  =  y  -  Yt 


X,  -  X, 

y.  -  y. 

Zf  -  z. 


(54) 


In  order  to  effectively  use  the  true  error,  its  magnitude,  |e,j,  is  compared 
with  the  rms  of  the  trace  of  Pg.  is  given  by 


P,  =  CPC^ 


(55) 


The  root  mean  square  of  the  trace  of  Pg  is  defined  as 

P-  =  (trPJ“  '  (oS  *  4  »  (56) 

where,  is  the  covariance  for  each  diagonal  position  element.  The  true  error  and 
rms  covariance  will  be  utilized  in  Section  3.3  and  following. 


Figure  7.  Performance  evaluation  with  no  feedback  of  a  Kalman  filter. 
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3.2  Truth  Model  Validation 


The  relative  position  part  of  the  truth  model  was  validated  by  running  a 
launch  trajectory  from  an  orbiting  vehicle  as  was  done  in  [3:56].  The  trajectory 
was  a  AV  of  60  m/sec  at  a  launch  angle  of  45°  from  a  near  circular  orbit.  The 
truth  model  orbit  was  given  an  eccentricity  of  0.00011.  Dunning  compared  the 
Clohessy-Wiltshire  solution  with  the  exact  solution  of  the  relative  equations  of 
motion.  Figure  8  shows  the  same  results  as  obtained  in  [3:56]. 


X 10*  TRUTHS  Relative  Trajectories 


Del  Y  [meters]  xlO^ 

Figure  8.  Comparison  of  the  Clohessy-Wiltshire  and  truth  model  solutions. 
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3.3  Filter  Model  Validation 


The  filter’s  general  behavior  is  determined  by  running  it  with  no  dynamics 
noise.  This  is  accomplished  by  providing  the  filter  with  =  0  and  analyzing  the 
behavior  of  the  true  error  and  the  rms  covariance.  With  =  0,  the  true  error, 
le,|,  and  rms  covariance,  P^,  were  plotted  and  are  shown  in  Figure  9. 

As  expected  the  true  error  diverges  while  the  rms  covariance  converges 
toward  zero  since  the  filter  has  total  confidence  in  its  dynamics.  This  confidence 
will  be  tempered  with  dynamics  noise  in  order  to  keep  the  filter’s  estimate  from 


Figure  9.  True  error  and  rms  covariance  for  on  and  off  with  no  dynamics 
noise. 
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reaching  perfection.  Notice,  however,  that  the  divergence  is  larger  in  the  case  with 
Jg  on.  This  also  is  expected  in  light  of  the  significant  non-linearities  and  coupling 
in  the  equations  of  motion  with  present. 

The  filter  was  tested  against  the  work  done  by  Capt  Johnston  [7:36-37], 
The  Jg  term  was  toggled  off  \n  order  to  directly  compare  the  results  with  his.  The 
Pq  and  from  Capt  Johnston’s  work  were  used  in  the  iterated,  extended  Kalman 
filter.  Results  were  similar  to  his.  After  minor  tuning,  the  iterated,  extended 
Kalman  filter  yielded  a  true  error  near  8  cm,  as  seen  in  Figure  10. 


J2  "ofr  Comparison 


Orbits 

Figure  10.  Jj  comparison  with  Capt  Johnston’s  work. 
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3.4  Filter  Tuning  Against  J2 

The  values  of  from  the  Jj  off  case  were  used  to  serve  as  a  starting  point 
to  tune  against  the  truth  model  with  Jg  on.  The  results  of  this  are  shown  in 
Figure  11.  As  can  be  seen,  more  tuning  is  needed  and  there  is  a  strong 
divergence  near  20  orbits.  After  numerous  attempts  to  adjust  by  keeping  the 
diagonal  elements  for  position  and  for  velocity  the  same,  the  filter  did  not  tune  out 
the  divergence  near  20  orbits.  This  was  investigated  further.  By  extending  to  86 
orbits  the  true  error  transient  at  20  orbits  appeared  to  be  the  only  transient.  The 


J2  on  Using  Tuned  J2  off  Qd 


I  Orbits 

Figure  11.  True  error  vs  rms  covariance  for  Jj  on  using  the  J2  offQ^. 
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behavior  of  the  filter  beyond  20  orbits  was  far  from  steady-state.  The  true  error 
was  diverging  as  can  be  seen  in  Figure  12.  The  transient  remamea  after  several 
attempts  to  achieve  steady-state  for  the  true  error  and  rms  covariance  by  adjusting 
the  position  elements  and  the  velocity  elements  the  same.  Since  the  filter  was 
running  in  single  precision,  the  possibility  of  a  numerical  precision  problem  was 
investigated. 
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3.4.1  Numerical  Precision  Investigation 


The  filter  was  programmed  to  operate  in  single  precision.  This  was  done 
to  simulate  the  on-board  software.  As  a  simple  check,  the  filter  was  run  in  double 
precision.  The  filter  performed  the  same  as  in  single  precision.  Figure  13  shows 
the  comparison  for  approximately  20  orbits.  Further  investigation  into  the 
dynamics  was  warranted. 
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3.4.2  Individual  Component  Investigation 


In  order  to  isolate  the  state  elements  that  were  causing  this  behavior,  the 
position  true  error  and  rms  covariance  were  plotted  for  each  axis.  Figure  14 
indicates  the  y  and  z  components  as  having  the  most  error  accounting  for  the 
transient  and  divergence.  The  Clohessy-Wiltshire  equations  uncouple  the  z 
position  and  velocity  components,  and  do  not  account  for  the  drift  in  the  relative 
orbits.  This  is  consistent  with  the  actual  relative  orbits  depicted  in  Figure  5. 
Following  this,  all  six  position  and  velocity  elements  were  then  adjusted 
individually. 

The  uncertainty  in  the  dynamics  is  attributed  to  an  uncertainty  in  the 
accelerations  due  to  Jg.  These  are  represented  by  velocity  elements  in  Q^.  To  get 
a  handle  on  these  elements,  an  acceleration  based  on  the  amount  of  drift  in  each 
direction  was  found.  From  these  accelerations,  the  elements  of  were  obtained. 

Appropriate  values  for  the  accelerations  were  obtained  from  the  distance  the 
vehicles  have  drifted  in  each  direction  per  orbit, 

(57) 

where  appropriate  values  for  were  obtained  from  Figure  5,  and  t  is  the  orbital 
period  of  approximately  6300  seconds.  Using  a  of  20  m  (200  m/10  orbits), 
^dhtt.y  approximately  1  X  10  m/s  .  Similarly,  for  of  5  m,  adnft./ X  10 
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Figure  14.  Component  breakout  of  true  error  and  rms  covariance. 
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m/s^.  The  x  component  shows  little  drift,  and  therefore  will  be  based  on  a 
^dntt.x  equal  to  the  range  measurement  error  assumed,  (.01  m).  This  equates  to  an 
^drin.x  of  5  X  10  ’°  m/s^. 

Using  these  accelerations,  better  estimates  of  the  velocity  elements  of 
were  found  by 


(  ^drin.x.y^^  0 


(58) 


Using  an  update  rate  of  300  seconds  yields  v^,  v^,  and  elements  of  2.25  X  10  '“* 
m^/s^,  9  X  lO"®  m^/s^,  and  5.6  X  10'^  m^/s^  respectively.  By  using  Eq  (57)  with  t  = 
At  a  better  estimate  of  the  position  elements  of  were  found.  Using  the  update 
rate  of  300  seconds  gives  the  x,  y,  and  z  elements  of  2.25  X  10  ®  m^  4.5  X  10^ 
m^,  and  1.1  X  10'^  respectively.  The  results  using  these  values  are  shown  in 
Figure  15. 

3.4.3  Utilizing  the  Kozai  Mean  Motion 

Even  with  these  values  for  Q^,  the  filter  still  was  not  tuning  out  the  transient. 
The  dynamics  model  in  the  estimator  was  augmented.  The  Clohessy-Wiltshire 
equations  are  in  the  desired  linear  form.  In  order  to  keep  the  equations  linear  the 
mean  motion  was  changed  to  the  Kozai  form. 
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J2  -on”  With  New  Qd 


Figure  15.  True  error  and  rms  covariance  with  a  better  estimate  of  the  velocity 
elements  of  Q^. 


The  Kozai  form  is  given  by  [4:369] 


3 

2 


(59) 


The  relative  orbits  displayed  in  Figure  16  show  the  effect  of  using 
This  approach  at  incorporating  Jj,  only  improves  the  prediction  of  the  drift  in  the 
y  direction.  It  doesn’t  couple  the  equations  of  motion.  Therefore,  the  y 
components  of  position  and  velocity  will  need  less  noise  than  anticipated  above. 
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Figure  16.  f?ea/ dynamics  and  Clohessy-Wiltshire  equations  using  (10  orbits). 
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With  this  in  mind  numerous  attempts  were  made  at  tuning  the  filter.  The  filter 
showed  results  that  had  some  sort  of  steady-state  behavior  as  seen  l:i  Figure  17. 


3.5  "Tuned"  Performance  Analysis 

The  filter  "tuned"  against  Jg  did  not  perform  as  well  as  against  the 
Newtonian  truth  model.  This  might  be  expected  since  the  filter  dynamics  are  a  l  st 
order  approximation  of  the  point  mass  orbit,  while  the  Jj  truth  model  is  highly  non- 
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Figure  17. 


J2  "on”  Performance 


"Tuned"  filter  performance. 


linear  and  coupled  with  a  perturbing  acceleration.  "Steady-state"  is  seen  in 
Figure  17  to  be  approximately  1  to  2  meters.  Although  its  performance  does  not 
compare  to  the  point  mass  analysis,  the  filter  still  remains  within  the  allowed 
tolerance  of  25  meters  for  the  86  orbits. 

3.6  Long  Term  Performance 

Although  Figure  17  shows  the  true  error  below  the  25  meter  accuracy 
requirement  for  as  long  as  86  orbits,  the  behavior  is  questionable  when  compared 
to  the  off  case.  The  final  analysis  was  to  determine  if  the  filter’s  dynamics  had 
sufficient  fidelity  to  model  the  effects  or  if  the  iterated,  extended  Kalman  filter 
was  not  providing  adequate  response.  In  order  to  assess  this,  both  the  Jg  on  and 
Jj  off  case  were  run  for  30  days. 

The  filter  remained  "steady-state"  for  quite  awhile  before  exhibiting  non¬ 
steady-state  behavior  in  the  true  error  for  the  Jg  on  case.  This  divergence  is 
greater  in  amplitude  than  the  previous  transient  at  20  orbits.  Figure  1 8  shows  the 
long  term  behavior  of  30  days.  The  true  error  diverged,  while  the  rms  covariance 
did  not  converge.  Figure  19  shows  the  off  case  run  for  30  days  as  well. 

Both  cases  were  run  using  sample  times  varying  from  112  seconds  to  300 
seconds  in  order  to  improve  the  filter’s  response.  No  significant  improvements 
resulted;  therefore,  the  300  second  sample  time  was  used  in  order  to  keep 
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[meters) 


computation  time  down. 

The  filter  showed  smooth  behavior  for  the  Jj  off  case  for  the  entire  30  days. 
This  indicates  the  iterated,  extended  Kalman  filter  has  adequate  ability  to  handle 
the  problem  if  the  dynamics  models  are  comparable.  Since  the  performance  in  the 
J2  on  case  is  poor  compared  to  the  Jj  off  case,  it  is  evidence  that  the  filter 
dynamics  do  not  have  adequate  ability  to  model  the  truth  dynamics. 
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IV.  Conclusions  and  Recommendations 


The  use  of  the  iterated,  extended  Kalman  filter  yielded  satisfactory  results 
for  the  Jj  off  case,  given  the  accuracy  requirements.  However,  the  analysis 
presented  here  indicates  the  performance  is  poor  for  the  filter  dynamics  model 
when  run  against  the  Jg  on  case.  This  Indicates  the  filter  dynamics  need  to  be 
modified  more  substantially  than  was  done.  The  Clohessy-Wiltshire  equations  had 
good  performance  against  the  Newtonian  point  mass  dynamics,  due  to  the  fact 
they  are  the  first  order  form  of  the  point  mass  dynamics.  Any  update  in  the  filter 
was  readily  propagated  forward  to  first  order  accuracy.  This  is  a  significant  factor 
of  the  filter’s  performance. 

However,  only  the  range  was  used  as  data,  and  therefore  only  partially 
modified  the  position  and  velocity  elements  of  the  state,  given  the  observation 
relationships,  Eqs  (41)  and  (44).  Against  the  Newtonian  point  mass  dynamics,  the 
velocity  elements  did  not  necessarily  need  to  be  updated  with  data  since  the 
dynamics  already  had  good  fidelity  of  the  truth  dynamics.  In  contrast,  against  the 
Jj  truth  dynamics,  the  fidelity  of  the  Clohessy-Wiltshire  equations  was  poor,  as 
seen  in  Figure  5  and  Figure  16.  This  is  primarily  due  to  the  uncoupling  of  the  z 
components  and  the  exclusion  of  Jj  effects. 

Although  the  Clohessy-Wiltshire  equations  are  valid  for  any  down  range 
relative  position,  they  are  valid  only  for  small  relative  positions  in  the  radial  and 
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out-of-plane  components  [13:80],  and  hence  only  small  periods  of  time.  The  filter 
uses  the  Clohessy-Wiltshire  trajectory  only  for  a  small  period  of  time,  the  sample 
period.  In  this  sample  period  the  state  changes  only  by  a  small  amount. 
Therefore,  the  dynamics  may  give  sufficient  performance  for  the  sample  period, 
but  the  state  not  being  fully  updated  may  prevent  adequate  performance  against 
J2  on.  The  significant  amount  of  tuning  by  adjusting  diagonal  elements  of  tried 
to  overcome  these  problems. 

Improved  performance  should  be  obtainable  by  either  improving  the  filter 
dynamics  or  by  including  range  rate  in  the  filter’s  observation  relationships.  The 
first  assumes  a  linear  form  of  relative  motion  equations  to  include  Jg  effects, 
eccentricity  and  coupling  of  the  equations  of  motion.  The  second  assumes  the 
satellite’s  ability  to  sense  range  rate.  However,  this  may  complicate  the  satellite’s 
hardware  design  too  much,  negating  the  idea  to  keep  the  satellites  as  simple  as 
possible. 

Further  work  should  include  a  better  filter  dynamics  model  such  as 
referenced  in  [6].  Alternately,  both  range  and  range  rate  observation  relationships 
should  be  used.  Having  established  the  new  dynamic  or  observation  relations,  a 
complete  Monte-Carlo  analysis  should  be  run.  Since  the  filter  is  required  to 
perform  for  several  satellites,  it  should  be  run  against  several  initial  conditions  to 
evaluate  its  robustness. 
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Appendix  A:  Haming  Subroutine 


SUBROUTINE  HAMING(NXT) 

Version  of  11/07/90 
Purpose 

Subroutine  for  integrating  a  system  of  first  order  differential 
equations.  It  is  a  fourth  order  predictor-corrector  algorithm 
which  means  it  carries  the  last  four  values  of  the  state  vector, 
and  extrapolates  these  values  to  obtain  a  predicted  next  value 
(the  prediction  step)  and  evaluates  the  equations  of  motion  at 
the  predicted  point,  and  then  corrects  the  extrapolated  point 
using  a  higher  order  polynomial  (the  correction  step). 

Input 

NXT  =  specifies  which  of  the  four  values  of  the  state  vector  is 
the  current  one.  NXT  is  updated  by  HAMING  automatically, 
but  must  be  set  to  ZERO  on  the  first  call. 

Call  Subroutines 

RHS(NXT)  =  evaluates  the  equations  of  motion 
External  Functions 
None 

Common  Blocks 

HAM  =  Memory  block  shared  by  the  main  driver  and  subroutine  RHS. 
The  common  block  contains: 

X  =  is  the  independent  variable  (often  time) 

Y(MAX,4)  =  the  state  vector  (4  copies),  with  NXT  pointing  to 

the  current  one,  the  limit  of  MAX  EOM  can  be  changed 
through  the  PARAMETER  in  main  driver,  sub  program 
RHS,  and  below. 

F(MAX,4)  =  are  the  EOM  evaluated  at  the  same  times  as  the  state 
vector  Y  ...  it  is  the  job  of  sub  program  RHS  to 
calculate  these. 

ERR(MAX)  =  is  an  estimate  of  the  one-step  integration  error 
N  =  is  the  number  of  ODES  ...  limit  is  MAX  unless  you  change 
the  PARAMETER  statement  in  main  driver,  sub  program 
RHS,  and  below. 

H  =  is  the  timestep  ...  one  call  to  HAMING  increments  X  by  H 
References 

"Numerical  Methods  for  Scientists  and  Engineers",  Richard  W.  Hamming; 


McGraw-  Hill,  2nd  Ed.,  1973;  pp.  361-408 


*  Donald  G.  M.  Anderson  -  Harvard  (1 972) 

*  Original  program  modified  by  Dr.  William  E.  Wiesel  and  Dr.  Rodney  Bain. 

*  Comments 

*  TOL  =  is  HAMING’s  startup  tolerance  ...  set  to  reasonable  value 

*  as  necessary  in  PARAMETER  statement. 

*  The  user  must  supply  a  main  driver,  and  the  subroutine  RHS(NXT) 

*  which  evaluates  the  equations  of  motion. 

IMPLICIT  REAL*8  (A-H,0-Z)  !  Global  double  precision 

PARAMETER  (ZERO=O.DO,  ONE=1.DO,  TWO=2.DO,  THREE=3.D0, 

1  FOU  R=4.  DO,  M AX=42,  TOL=  1 .  D- 1 2) 

COMMON  /HAM/  X,Y(MAX.4).F(MAX,4).ERR(MAX).N.H 

*  Check  if  this  is  the  first  call  ...  HAMING  (like  all  predictor- 

*  correctors)  needs  ’previous’  values 

IF(NXT)  190,10,200 

*  It  is  a  forward  Picard  iteration  (slow  and  expensive)  to  step 

*  backwards  in  time  three  steps  to  get  the  4  previous  points.  A 

*  successful  startup  returns  NXT=1,  and  time  has  not  been 

*  incremented.  If  startup  fails,  NXT  will  be  returned  as  ZERO. 

10  XO=X 
HH=H/TWO 
CALL  RHS(1) 

DO  40  L=2,4 
X=X+HH 
DO  20  1=1, N 

20  Y(I,L)=Y(I,L-1)+HH*F(I,L-1) 

CALL  RHS(L) 

X=X+HH 
DO  30  1=1, N 

30  Y(I,L)=Y(I,L-1)+H*F(I,L) 

40  CALL  RHS(L) 

JSW=-10 
50  ISW=1 

DO  120  1=1, N 

HH=Y(I,1)+H*(9.D0*F(I,1)+19.D0*F(I.2)-5.D0*F(I,3) 


50 


1  +F(I,4))/24.D0 

IF(  DABS(HH-Y{I,2))  .LT.  TOL)  GOTO  70 
ISW=0 

70  Y(I,2)=HH 

HH=Y(I.1  )+H*(F(l.1  )+FOUR*F(l.2)+F(l.3))rrHREE 
IF(  DABS(HH-Y(I.3))  .LT.  TOL)  GOTO  90 
ISW=0 

90  Y(I,3)=HH 

HH=Y(I,1)+H*(THREE*F(I.1)+9.D0*F(1,2)+9.D0*F(I,3) 

1  +THREE*F(I.4))/8.D0 
IF(  DABS(HH-Y(I.4))  .LT.  TOL)  GOTO  110 
ISW=0 

110  Y(I.4)=HH 
120  CONTINUE 
X=XO 

DO  130  L=2,4 
X=X+H 

130  CALLRHS(L) 

IF(ISW)  140,140,150 
140  JSW=JSW+1 
IF(JSW)  50,280,280 
150  X=XO 
ISW=1 
JSW=1 
DO  160  1=1  ,N 
160  ERR(l)=ZERO 
NXT=1 
GOTO  280 

*  A  call  to  NAMING  with  NXT=-NXT,  after  a  successful  startup. 

*  will  turn  off  the  second  evaluation  of  the  equations  of  motion 

*  following  the  corrector  step.  In  systems  where  the  equations  of 

*  motion  are  very  expensive,  this  can  halve  your  run  time. 

190  JSW=2 

NXT=IABS{NXT) 


*  This  is  the  predictor-corrector  algorithm  ...  first  the  indices 

*  are  premuted. 

200  X=X+H 
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NP1=M0D(NXT.4)+1 
GOTO  (21 0.230), ISW 
210  GOTO  (270,270,270,220).NXT 
220  ISW=2 

230  NM2=MOD(NP1,4)+1 
NM1=MOD(NM2,4)+1 
NPO=MOD(NM1,4)+1 

* ...  then  the  predictor  part  is  run  to  find  an  extrapolated  value 

*  of  the  state  vector  at  the  new  time  ... 

DO  240  1=1, N 

F(I,NM2)=Y{I.NP1  )+FOUR*H*(TWO*F(l,NPO)-F{l,NM1 ) 

1  +TWO*F(l,NM2))/THREE 

240  Y(I.NP1  )=F(I,NM2)-0.92561 9835D0*ERR(I) 

*  The  equations  of  motion  are  evaluated  at  the  extrapolated  value 

*  of  the  state  vector ... 

CALL  RHS(NPI) 

*  and  the  corrector  algorithm  is  used  to  add  this  new  information 

*  and  obtain  a  better  value  of  the  new  state  vector ... 

DO  250  1=1  ,N 

Y(I,NP1  )=(9.DO*Y(I.NPO)-Y(I.NM2)+THREE*H*(F{I,NP1 ) 

1  +TWO*F(I,NPO)-F(I,NM1)))/8.DO 

ERR(I)=F{I,NM2)-Y(I,NP1) 

250  Y(I,NP1  )=Y(I,NP1  )+0.0743801 653D0*ERR{I) 

GOTO  (260,270),JSW 

*  Finally,  the  equations  of  motion  are  re-evaluated  at  the  better 

*  value  of  the  state  vector ...  this  can  be  suppressed. 

260  CALLRHS(NPI) 

270  NXT=NP1 

280  RETURN 
END 


52 


Bibliography 


1.  Bate,  Roger,  Donald  D.  Mueller,  and  Jerry  E.  White.  Fundamentals  of 
Astrodynamics.  New  York:  Dover  Publications,  1971. 

2.  Boden,  Daryl  G.  and  Bruce  A.  Conway.  "A  Comparison  of  Nonlinear  Filters 
for  Orbit  Determination",  Astrodynamics  Conference,  AIAA  Paper  86-2055. 
Williamsburg,  VA:  American  Institute  of  Aeronautics  and  Astronautics,  August 
1986. 

3.  Dunning,  Robert  Scott.  "The  Orbital  Mechanics  of  Flight  Mechanics".  NASA 
SP-325.  National  Aeronautics  and  Space  Administration,  1 973. 

4.  Escobal,  Pedro  R.  Methods  of  Orbit  Determination.  John  Wiley  and  Sons, 
Inc.,  1965. 

5.  Filer,  Captain  Sherrie  Norton.  Investigation  of  the  Observability  of  a  Satellite 
Cluster  in  a  Near  Circular  Orbit.  MS  Thesis,  AFIT/GA/ENY/89D*2,  School  of 
Engineering,  Air  Force  Institute  of  Technology  (AU),  Wright-Patterson  AFB 
OH,  December  1989. 

6.  Hujsak,  R.  S.  "A  Non-linear  Dynamical  Model  of  Relative  Motion  for  the 
Orbiting  Debris  Problem."  Astrodynamics  1989:  Proceedings  of  the 
AAS/AIAA  Astrodynamics  Conference  (AAS/AIAA  Paper  89-397). 
American  Astronautical  Society,  San  Diego  CA,  1989. 

7.  Johnston,  Captain  Stephen  C.  Autonomous  Navigation  of  a  Satellite  Cluster. 
MS  Thesis,  AFIT/GA/ENY/90D-9,  School  of  Engineering,  Air  Force  Institute 
of  Technology  (AU),  Wright-Patterson  AFB  OH,  December  1990. 

8.  Maybeck,  Peter  S.  Stochastic  Models,  Estimation,  and  Control  Volume  1. 
New  York:  Academic  Press,  1 979. 

9.  Maybeck,  Peter  S.  Stochastic  Models,  Estimation,  and  Control  Volume  2. 
New  York:  Academic  Press,  1 982. 

1 0.  Ward,  Captain  Michael  L.  P.  Estimated  Satellite  Cluster  Elements  in  Near 
Circular  Orbit  MS  Thesis,  AFIT/GA/AA/88D-13,  School  of  Engineering,  Air 
Force  Institute  of  Technology  (AU),  Wright-Patterson  AFB  OH,  December 
1988. 


53 


1 1 .  Wiesel,  William  E.  Class  text  distributed  iii  MECH  731 ,  Modern  Methods  of 
Orbit  Determination.  School  of  Engineering,  Air  Force  Institute  of  Technology 
(AU),  Wright-Patterson  AFB  OH,  April  1990. 

12.  Wiesel,  William  E.  Class  text  distributed  in  MECH  636,  Advanced 
Astrodynamics.  School  of  Engineering,  Air  Force  Institute  of  Technology  (AU), 
Wright-Patterson  AFB  OH,  January  1 991 . 

13.  Wiesel,  William  E.  Spaceflight  Dynamics.  New  York:  McGraw-Hill,  1989. 


54 


REPORT  DOCUMENTATION  PAGE 


form  Approved 
0MB  No  0704  0188 


P'jOJic 'eDorT'"9  Ourae^  ;:iier’ion  :f  nfcrmatrcr  s  i.e'vae  ' -ouf  oe^' 'e^oorse.  ncuamg  th#  r.me  *C'' '■evewing  rnstr^^aicn?  r^•- j  ,our':in 

gathering  and  maintaining  the  data  neeoeo.  ana  comoieting  ann  '•eviev^.ng  tne  rr-iferr-on  ot  ntormation  Sena  comments  regarding  this  'xi^aen  est  mate  :r  jn*  ;tner  iboect  'nn 
coileaion  of  ntormation  nciuamg  suggestions  ’cr  reducing  tn.>  ouroen  ’o  rv.jsningtor'  '^eaaQuar»ers  sen^Kes.  Directorate  ’or  ntormation  Ooeraricns  ana  =**»oc'^s.  f ^  e'^e'-son 
Davis  Highway.  Suite '204  Arlington,  .a  22202-4302  and  to  me  Ot*ie  j' v^anagement  and  Budget  ®aoerwcr<  Recuctior  P'oie.t  lC7C4*0'%8).  Washington  ;c  2012  3 


1.  AGENCY  JSE  ONLY  (Leave  Olank) 


2.  REPORT  DATE 

3  REPORT  TYPE  AND  OATES  COVERED 

December  1901 

Master's  Thesis 

4.  TITLE  AND  SUBTITLE 

S.  FUNDING  NUMBERS 

MA'/IGATION  OF  A  SATELLITE  CLL’STER  WITH  REALISTIC  DYNAMICS 

6.  AUTHOR(S) 

J.  Timothy  Middendorf 

7.  PERFORMING  ORGANIZATION  NAME($)  AND  ADORESS(E$) 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

Air  Force  Institute  of  Technology 

WPAFB,  OH  45433-6583 

AFIT/GA/F.NY7  91D-3 

9.  SPONSORING  MONITORING  AGENCY  NAME(S)  AND  AODRESS(ES) 

10.  SPONSORING  MONITORING 

AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES 

12a.  DISTRIBUTION  AVAILABILITv  STATEMENT 

12b.  DISTRIBUTION  CODE 

Approved  tor  public  release;  distribution  unlimited 

13.  ABSTRACT  {Maximum  200  wordi) 

Previous  work  in  the  area  of  estimation  of  relative  positions  within  a  satellite 
cluster  shoved  favorable  results.  However,  the  work  was  done  using  point  mass 
orbits  in  the  truth  model.  This  thesis  investigates  the  estimation  of  relative 
satelli-e  positions  operating  in  near  circular  orbitsinc luding  the  J2  term  in. 
Earth's  geopotential.  The  iterated,  extended  KAlman  filter  is  used  as  the  on-board 
estimator  in  order  to  gain  better  performance  in  the  face  of  the  non-linearities. 
The  dynamiv-S  in  the  estimator  are  based  on  the  Clohessy-Wi Itshire  equations  for 
relative  orbital  motion.  Inputs  to  the  host  estimator  are  range  measurement.'  from 
each  satellite  in  the  cluster.  The  relative  position  to  the  host  satellite  is 
investigated.  A  comparison  of  the  true  error  and  the  root  mean  squar  covariance 
was  performed . . 


14.  SUBJE..T  TERMS 

Satellite  Cluster  K 

man  Filter  Clohessy-Wi Itshire 

15.  NUMBER  OF  PAGES 

67 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFIC*";CN 

18.  SECURITY  CLASSIFICATION 

19.  SECURITY  CLASSIFICATION 

20.  LIMITATION  OF  ABSTRACT 

OF  REPORT 

OF  THIS  PAGE 

OF  ABSTRAa 

Unclassified 

Unclassified 

Unclassified 

UL 

NSN  75AO-01  280-5500 


Standard  -orm  298  °e~i  '  ?9 

P'WCriOW  Ov  A%S  •.*  . 

298-’C2 


